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ABSTRACT 

Intense structure formation and reionization occur at high redshift, yet there 
is currently little observational information about this very important epoch. 
Observations of gravitational waves from massive black hole (MBH) mergers 
can provide us with important clues about the formation of structures in the 
early universe. Past efforts have been limited to calculating merger rates using 
different models in which many assumptions are made about the specific values of 
physical parameters of the mergers, resulting in merger rate estimates that span 
a very wide range (0.1 — 10 4 mergers/year). Here we develop a semi-analytical, 
phenomenological model of MBH mergers that includes plausible combinations 
of several physical parameters, which we then turn around to determine how 
well observations with the Laser Interferometer Space Antenna (LISA) will be 
able to enhance our understanding of the universe during the critical z ~ 5 — 30 
structure formation era. We do this by generating synthetic LISA observable data 
(total BH mass, BH mass ratio, redshift, merger rates), which are then analyzed 
using a Markov Chain Monte Carlo method. This allows us to constrain the 
physical parameters of the mergers. We find that our methodology works well 
at estimating merger parameters, consistently giving results within 1-cr of the 
input parameter values. We also discover that the number of merger events 
is a key discriminant among models. This helps our method be robust against 
observational uncertainties. Our approach, which at this stage constitutes a proof 
of principle, can be readily extended to physical models and to more general 
problems in cosmology and gravitational wave astrophysics. 
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1. Introduction 


The non-linear universe is a new frontier in cosmology, especially at relatively high 
redshifts (z ~ 5 — 30) because that is when baryonic structure formation begins (e.g., Barkana 
& Loeb 2001). The first objects to form, dark matter halos, undergo hierarchical mergers 
and assemble into larger structures (Rees & Ostriker 1977; White & Rees 1978; Peebles 
1993). The first galaxies formed inside these halos; however, these primordial galaxies (with 
M ~ 1O 8 M 0 ) are very dim and hence not readily observed. Thus we are (quite literally) 
left in the dark regarding this crucial era in the history of the universe. Therefore, new 
techniques must be developed that will allow us to probe the physical processes that drive 
the formation of structures in the early universe. One such technique involves studying 
the mergers of massive black holes (MBHs), which trace the mergers of their parent halos. 
MBH mergers can be observed at high redshifts even if their host galaxies are not bright 
enough to be detected (e.g., Volonteri et al. 2003; Wyithe & Loeb 2003; Sesana et al. 2004; 
Hughes & Menou 2005). Consequently, observations of the gravitational waves emitted by 
MBH mergers at high- £ can further our understanding of structure formation in the early 
universe. 

Due to their importance in the study of structure formation, many research groups 
have calculated merger rates of MBHs (e.g., Haehnelt 2003; Enoki et al. 2004; Islam et al. 
2004; Sesana et al. 2004, 2005; Rhook & Wyithe 2005). However, the results span a very 
broad range (0.1 — 10 4 mergers/year). The large discrepancies that arise in the assessment 
of merger rates result from the choice of different physical parameters in the models used to 
perform the calculations. Preliminary results suggest that the minimum mass for a halo to 
host a black hole is of special importance in the determination of merger rates. 

Here we develop a semi-analytical, phenomenological model that includes plausible com- 
binations of several physical parameters involved in MBH mergers. As a proof of principle 
we use statistical methods to generate synthetic LISA observable data (total BH mass, BH 
mass ratio, redshift, merger rates), which are then run through a Markov Chain Monte Carlo 
(MCMC) algorithm to constrain the physical parameters associated with the mergers. Be- 
cause LISA distributions of observable parameters are expected to be broad, our method is 
relatively robust against uncertainties in the observations. Indeed, our results suggest that 
even if there are measurement biases, the total number of L/SM-detected events will be a 
strong discriminant among models. Our method is also general enough that it can be used 
with a number of different observational data constraints. 

In § 2 we establish a framework for understanding the dynamics of halo mergers using the 
standard Press-Schechter and Extended Press-Schechter approaches. In § 3 we discuss the 
relationship between MBHs and their host halos. This is the source of most of the uncertainty 
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in merger models. We present our simulations and results in § 4, then summarize in § 5. 


2. Dynamics of Halo Mergers 

The theory of hierarchical structure formation is based on the growth of linear pertur- 
bations in an initial cold dark matter (CDM) density field. Specifically, the Press-Schechter 
formalism (Press & Schechter 1974) gives the number density of dark matter halos as a func- 
tion of mass and redshift. We use the Sheth & Torrnen (1999) modification of the original 
Press-Schechter theory in our calculations, as it provides a better fit with numerical simu- 
lations. A nicely detailed summary of these standard procedures is given by Mo & White 
( 2002 ). 

The Press-Schechter formalism does not, however, allow us to follow the merger history 
of any particular halo. For this we use the Extended Press-Schechter formalism of Lacey 
& Cole (1993). This formalism has been criticized because of mathematical inconsistencies 
that arise in the calculations of merger trees (Somerville & Kolatt 1999; Benson et al. 2005). 
However, it is important to note that we do not construct full merger trees. Instead, we 
recalculate the Press-Schechter halo distribution at each redshift interval. This approach 
eliminates cumulative errors, since the inconsistencies disappear for small redshift intervals. 

Throughout our calculations we use cosmological parameters determined by the WMAP 
3rd- year results (Spergel et al. 2007). Although newer results are available, our proof of 
principle method does not need such high precision to work properly. 


3. The Relationship Between Massive Black Holes and their Host Halos 

The largest uncertainty in merger rate calculations involves the relationship between 
MBHs and their host dark matter halos. Here we describe some of these issues. 


3.1. Minimum Halo Mass and the Black Hole Occupation Fraction 

Suppose no black holes can form in dark matter halos that have less than some minimum 
mass M min . We motivate this with the standard assumption that black holes need to form 
from baryons. That is because baryons can cool, allowing for the formation of stars which 
can then evolve into black holes. If the baryons cannot be kept within the halo, then this 
does not happen and no black holes can form. The minimum mass of a halo that can host a 
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black hole therefore depends on the depth of the potential well compared to the temperature 
of the baryons. 

There are various possible mechanisms for the formation of massive black holes (Rees 
1984). The most common scenario is the collapse of Population 111 stars (Abel et al. 2002; 
Bromm et al. 2002; Heger et al. 2003). These are very massive (~ 100M e ) metal- free stars 
which form at high redshifts in halos with masses 10 5 — 1O 6 M 0 . Thus halos that contain 
Population 111 stars are likely to be populated with a MBff after the star dies. Here we 
assume for simplicity that all halos with mass above a certain M min are already occupied 
with black holes by z = 20, and that halos with mass below M min have a decreasing non-zero 
probability of being occupied by a black hole, P occ = (M/M min ) p , with p > 0. 


3.2. Black Hole-Halo Mass Correlation 

Even if we do have a halo that can form a black hole, there is uncertainty about the 
relation of the mass of the central black hole to the mass of the dark matter halo. Though 
a link has been shown between the central velocity dispersion of stars and the mass of a 
supermassive black hole (see Ferrarese & Ford 2005 for a review), there has been less work 
establishing a link between dark matter halo mass and black hole mass. However, Ferrarese 
(2002) found the following relation in the local universe: 

This can be generalized to include a redshift dependence: 

M b „ = 10 7 M s (1 + 2 )”, (2) 

where scaling arguments suggest n — 5/2 (e.g., Rhook & Wyithe 2005). 


3.3. Other issues not treated in the current model 

There are various other issues which are important in a full physical model of black 
hole mergers. However, these are not essential in our proof of principle, since our goal here 
is to demonstrate the method. Thus we make certain assumptions about them without 
parameterizing. Two of these issues, which will be treated in more detail in later work, are 
dynamical friction and gravitational recoil. 
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In order for a black hole merger to occur, there must first be a merger between two halos 
with black holes at their centers. Dynamical friction plays an important role in this process, 
as its drag brings the two halos together. Once the halos have merged, dynamical friction 
causes the central black holes to get closer and form a binary. Note, however, that dynamical 
friction between halos becomes inefficient for very large halo mass ratios. Therefore, in our 
work we assume that for any two halos with masses M\ and M 2 , where M\ > M 2 , there is 
a mass ratio cutoff Mi/M 2 = 50 above which no mergers happen. Additionally, if a merger 
takes too long, there is the probability of a third object coming in and disrupting the binary, 
thus ejecting one or more black holes in the process (Hoffman & Loeb 2007). Thus we 
assume that for M\jM 2 < 50 mergers happen immediately, thus not allowing for three-body 
interactions. 

After a black hole binary hardens, gravitational wave emission will make the black holes 
plunge and merge. However, if the black holes have unequal masses or spins there will be 
a gravitational recoil effect (e.g., Baker et al. 2006; Campanelli et al. 2007; Herrmann et 
ah 2007; Koppitz et ah 2007; Briigmann et ah 2008) that can eject the merger remnant 
out of the halo if the kick velocity is high enough (Madau & Quataert 2004). In our work 
we assume that if the kick speed of a merger is greater than the escape speed of the host 
halo (r’kick > v esc ), then there is a high probability of ejection of the merger remnant. For 
simplicity, we consider only non-spinning black holes. Black holes with spin will be treated 
in a future work. 


4. Simulations and Results 

We want to answer the question of what would LISA observations tell us about the 
process of structure formation in the early universe. To do this, we generate synthetic LISA 
observable data (total BH mass, BH mass ratio, redshift, merger rates) based on a simple 
model of MBH mergers that involves various assumptions about a handful of merger param- 
eters. We choose input values for these parameters within a range of plausible, physically 
realistic extremes. Using these inputs, we calculate the number density of halos using the 
procedures described in § 2. We then sample observables from the probability distribution 
thus generated by using the rejection method (Gentle 2003). The generated number of events 
is determined by the observation time. For our simulations the assumed observation time is 
three years. 

Our statistical goal is to find and explore regions of high likelihood in parameter space. 
The synthetic data is analyzed to determine how well we can recover the values of the input 
parameters used in the data generation process. The Markov Chain Monte Carlo (MCMC) 
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method gives us a computationally robust yet inexpensive way of doing this. A useful 
description of the general MCMC algorithm is given by Verde et ah (2003). 

We are currently working with a four-parameter model of merger events. These param- 
eters, and their explored ranges, are: 

(i) Minimum halo mass (M min ). This is the mass above which all halos are occu- 
pied with black holes at z = 20. The explored range of this parameter is 5.0 < 
log lo (M min /M 0 ) < 12.0. 

(ii) Power law index ( p ). Indicates the probability of occupation of a halo with mass 
M < M min at z = 20 as P occ = (M / M min ) p . for 0.1 < p < 4.0. 

(iii) Redshift dependence (n). Indicates the redshift dependence of the BH-halo mass 
relationship, and it is given by M B h,o = [M BH (M halo )](l + z) n , for 0.0 < n < 5.0. 

(iv) BH mass spread (cr). This is the spread in the BH mass as a function of halo mass, 
defined as P(log 10 (M BH /M o )) oc exp[(log 10 (M BH /M o ) - log lo (M BH , o /M 0 )) 2 /2(T 2 ], for 
0.1 < a < 1.0. 

We generated synthetic data for five combinations of input parameter values. The top 
portion of Table 1 shows the input values and the results of the simulations, and Figure 1 
illustrates a representative sample. Additionally, we simulated the effects of observational 
uncertainty by adding errors to the synthetic LISA data. This is done by scrambling the 
values of the total BH mass, BH mass ratio and redshift independently of each other. That 
is, for some parameter x with unaltered value ay, we created a biased data set by altering 
the parameter value to (1.0 — 1A)xq, and an unbiased data set by altering the parameter 
value to (0.8 — 1.2)x 0 . The bottom portion of Table 1 shows the results of these tests, and 
Figure 2 shows a representative sample. 

Table 1 (top section) shows the parameter values at the point of maximum likelihood 
obtained from our MCMC procedure in each simulation. The error bars represent the 68% 
confidence level, or a l-cr deviation. A quick comparison between the results and the input 
parameter values reveals that for Runs 1, 2, 3 and 5 our points of maximum likelihood fall 
within l-cr of their input value, while Run 4 has one resultant parameter value lie just outside 
the l-cr margin. Of the 20 parameter ranges explored (five simulations with four parameters 
each), only one is outside the 68% confidence level. This suggests a slightly non- Gaussian 
distribution. 

Figure 1 illustrates the results from one of these simulations (Run 1). Each plot 
shows the relationship between a pair of parameters (left column: log lo (M min /M 0 ) — p, 
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log lo (M min /M 0 ) — n, log lo (M min /M 0 ) — a; right column: p — n, p — a, n — a). Notice the 
nicely constrained correlations, particularly in the log lo (M min /M 0 ) — p plot. This is because 
the number of merger events is critical in determining the shape of the confidence regions. 

The bottom part of Table 1 (Runs lb, lu, 2b, 2u) shows the parameter values at the 
point of maximum likelihood obtained from our MCMC procedure on each simulation done 
with data altered in a biased or unbiased way. From the table data it can quickly be assessed 
that simulations carried out with altered unbiased data give us results as accurate as those 
for which the data were not altered, i.e., all resultant parameter values lie within 1-cr of their 
input values. However, for altered biased data we can see that the results lie outside the 1-cr 
region. 

Figure 2 shows the results of a representative altered-data simulation. All three plots 
show the log lo (M min /M 0 ) — p correlation; the top panel resulting from altered unbiased 
data (Run 2u) and the bottom panel resulting from altered biased data (Run 2b). The 
middle panel (Run 2) shows the results from the control simulation, for which the data was 
not altered in any way. Notice that the results from altered unbiased data keep the point 
of maximum likelihood within the 68% confidence level, while the results from the altered 
biased data do not. Additionally, altered unbiased data results in a wider 1-cr region than 
unaltered data. Note however that we obtain a fairly reasonable parameter estimation even 
from altered unbiased data, and biased data as well but to a lesser extent, because in our 
simple scenario the number of merger events can discriminate between models with different 
parameter values. 


5. Summary and Discussion 

Our method works well at estimating merger parameters even when there are observa- 
tional uncertainties involved. Our simulations run with altered unbiased data give us results 
of comparable accuracy to the results obtained from simulations where the data was unal- 
tered, albeit with slightly less tight constraints. Results from simulations run with altered 
biased data are expectedly less accurate, with some parameter values estimated outside of 
the 1-cr margin from the input values. The general shape of the confidence regions, however, 
remains consistent throughout our results, even when altered biased data was used. This 
is because the number of merger events determines the shape of the confidence region, thus 
allowing us to distinguish models with different parameter values. 

We have thus established a flexible framework which works well at estimating MBH 
merger parameters in our phenomenological model. In future work we will explore more 



physically driven models with additional merger and cosmological parameters. One espe- 
cially important parameter is the redshift of reionization. Reionization can increase the 
temperature of the baryons at the centers of dark matter halos, effectively resulting in an 
increase in the minimum mass for a halo to form a black hole. Another essential model com- 
ponent is the LISA detection sensitivity, since not every merger that happens in the universe 
is detectable. This can also lead to degeneracies between models, i.e., very different models 
resulting in similar numbers of detected merger events (Sesana et al. 2007). Furthermore, 
the statistical method we use is sufficiently general that in the future it will be possible to 
incorporate additional sources of parameter constraints, such as currently available observa- 
tional data of high redshift galaxies and future observations carried out with the Atacama 
Large Millimeter Array (ALMA) and the James Webb Space Telescope ( JWST ) once they 
become operational. 
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EAM was supported by the NASA/GSFC Cooperative Education Program and MCM was 
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Table 1. Simulations 


Input Parameters Max Likelihood 


Run 

N ev 

m 

P 

n 

<7 

m 



P 



n 



o 


1 

169 

8.5 

1.5 

0.5 

0.2 

8.62 

± 

0.36 

1.42 

± 

0.45 

0.55 

± 

0.20 

0.19 

± 

0.10 

2 

245 

9.3 

1.0 

1.5 

0.4 

9.43 

± 

0.36 

0.94 

± 

0.11 

1.68 

± 

0.27 

0.39 

± 

0.08 

3 

219 

9.3 

1.0 

1.5 

0.2 

9.13 

± 

0.28 

1.07 

± 

0.12 

1.59 

± 

0.23 

0.22 

± 

0.08 

4 

16 

9.3 

1.5 

1.5 

0.4 

9.08 

± 

0.46 

1.75 

± 

1.36 

1.45 

± 

0.54 

0.18 

± 

0.21 

5 

155 

9.3 

1.0 

0.5 

0.4 

9.19 

± 

0.46 

1.07 

± 

0.21 

0.68 

± 

0.28 

0.39 

± 

0.09 

lb 

169 

8.5 

1.5 

0.5 


8.54 

± 

0.24 

1.57 

± 

0.35 

0.99 

± 

0.17 

0.11 

± 

0.07 

lu 

169 

8.5 

1.5 

0.5 


8.63 

± 

0.36 

1.41 

± 

0.41 

0.52 

± 

0.20 

0.22 

± 

0.10 

2b 

245 

9.3 

1.0 

1.5 


8.93 

± 

0.21 

1.16 

± 

0.12 

1.85 

± 

0.20 

0.35 

± 

0.07 

2u 

245 

9.3 

1.0 

1.5 


9.51 

± 

0.38 

0.92 

± 

0.12 

1.74 

± 

0.28 

0.39 

± 

0.08 


Note. — Col. (1): Simulation number. Col. (2): Number of generated merger events 
over a three year period. Cols. (3, 7): Minimum halo mass, m = log lo (M min /M 0 ). Cols. 
(4, 8): Power law index, p. Cols. (5, 9): Redshift dependence, n. Cols. (6, 10): BH 
mass spread, a. Columns 3 — 6 are the input parameters for the simulations. Columns 
7 — 10 show the maximum likelihood results for each simulation (i.e., the best estimate 
parameters). Error bars represent the 68% confidence level. Simulations labeled ’b’ and ’u’ 
represent biased and unbiased data, respectively. Results show our method works well at 
estimating merger parameters. 
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i A 95% 
a □ 90% 



Fig. 1. — Confidence regions in parameter estimation for 169 merger events over a 3-year 
observation period. The actual values of the parameters are log 10 (M min /M o ) = 8.5, p = 1.5, 
n = 0.5, a = 0.2 (Run 1). Notice the nicely constrained correlation on the first panel 
(top- left). 
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Fig. 2. — Simulating observational errors. Top panel: altered unbiased data (Run 2u). 
Middle panel: no alterations (Run 2). Bottom panel: altered biased data (Run 2b). Actual 
values of the parameters shown here are log 10 (M min /MQ) = 9.3 and p = 1.0, with 245 total 
merger events over a 3-year observation run. Note that the parameter estimation is not 
severely affected by either biased or unbiased errors in the data. 


